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Electronic Structure of Multilayer Graphene 

Hongki MlN*) and Allan H. MacDonald 
Department of Physics, University of Texas at Austin, Austin, Texas 78712, USA 

We study the electronic structure of multilayer graphene using a 7r-orbital continuum 
model with nearest-neighbor intralayer and interlayer tunneling. Using degenerate state 
perturbation theory, we show that the low-energy electronic structure of arbitrarily stacked 
graphene multilayers consists of chiral pseudospin doublets with a conserved chirality sum. 

§1. Introduction 

The recent explosiorllHl} of research on the electronic properties of single layer 
and stacked multilayer graphene sheets has been driven by advances in material 
preparation methods JH 1 '® by the unusuaHMU'JZl electronic properties of these mate- 
rials including unusual quantum Hall effects J^'FH and by hopes that these elegantly 
tunable systems might be useful electronic materials. 

In this paperj!2f we study the electronic structure of arbitrarily stacked multi- 
layer graphene using a 7r-orbital continuum model with only near-neighbor interac- 
tions, analyzing its low-energy spectrum using degenerate state perturbation theory. 
Here we focus solely on aligned multilayer graphene without rotational stacking 
feultsED Interestingly, we find that the low-energy effective theory of multilayer 
graphene is always described by a set of chiral pseudospin doublets with a conserved 
chirality sum. We discuss implications of this finding for the quantum Hall effect in 
multilayer graphene. 

§2. 7r-orbital continuum model 

We consider the 7r-orbital continuum model for iV-layer graphene Hamiltonian 
which describes bands near the hexagonal corners of the triangular lattice Brillouin 
zone, the K and K' points: 

p 

where <Jp = (ci iQ)P , cx s p jP , ■ ■ ■ , CN, a ,pi c iV,/3, P ) and Q iMjP is an electron annihilation 
operator for layer I = 1, ••■ ,N, sublattice \x = a, (3 and momentum p measured 
from K or K' point. 

The simplest model for a multilayer graphene system allows only nearest-neighbor 
intralayer hopping t and the nearest-neighbor interlayer hopping t±. The in-plane 
Fermi velocity v is related with t by ~ = ^t, where a = 2.46 A is a lattice constant 
of monolayer graphene. Although this model is not fully realistic, some aspects of 
the electronic structure can be understood by fully analyzing the properties of this 
simplified model first and then considering corrections. 
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Fig. 1. (Color online) (a) Energetically favored stacking arrangements for graphene sheets. The 
honeycomb lattice of a single sheet has two triangular sublattices, labeled by a and (3. Given 
a starting graphene sheet, the honeycomb lattice for the next layer is usually positioned by 
displacing either a or f3 sublattice carbon atoms along a honeycomb edge. There are therefore 
in three distinct two-dimensional (2D) sheets, labeled by A, B, and C. Representative a and (3 
sublattice positions in A, B, and C layers are identified in this illustration. It is also possible 
to transform between layer types by rotating by ±60° about a carbon atom on one of the two 
sublattices. (b) Each added layer cycles around this stacking triangle in either the right-handed 
or the left-handed sense. Reversals of the sense of this rotation tend to increase the number 
of low-energy pseudospin doublets Nd- In graphite, Bernal (AB) stacking corresponds to a 
reversal at every step and orthorhombic (ABC) stacking corresponds to no reversals. 



2.1. Stacking diagrams 

When one graphene layer is placed on another, it is energetically favorabldlD for 
the atoms of either a or (3 sublattices to be displaced along the honeycomb edges, 
as illustrated in Fig. [TJ This stacking rule implies the three distinct but equivalent 
projections (labeled A, B, and C) of the three-dimensional structure's honeycomb- 
lattice layers onto the x-y plane and 2 N ~ 2 distinct iV-layer stack sequences. When 
a B layer is placed on an A layer, a C layer on a B layer, or an A layer on a C layer, 
the a sites of the upper layer are above the (5 sites of the lower layer and therefore 
linked by the nearest interlayer neighbor 7r-orbital hopping amplitude t±. For the 
corresponding anticyclic stacking choices (A on B, B on C, or C on A), it is the (5 
sites of the upper layer and the a sites of the lower layer that are linked. All distinct 
N = 3, N = 4, and N = 5 layer stacks are illustrated in Fig. [2] in which we have 
arbitrarily labeled the first two layers starting from the bottom as A and B. 

2.2. Energy band structure 
2.2.1. Preliminaries 

Before analyzing energy spectrum of multilayer graphene, let us consider the 
Hamiltonian of a one-band tight-binding model for a chain of length N with near- 
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Fig. 2. (Color online) Stacking sequences and linkage diagrams for iV = 3, 4, 5 layer stacks. The 
low-energy band and Landau level structures of a graphene stacks with nearest-neighbor hopping 
are readily read off these diagrams as explained in the text. Shaded ovals link a and (5 nearest 
interlayer neighbors. 



neighbor hopping parameter t±: 





( ° 


t± 












t± 





t± 







H = 





t± 





*1 • 












t± 








V 



(2-2) 



/ 



This Hamiltonian is important for analyzing the role of interlayer hopping as we 
explain below. 

Let a = (ai, ajv) be an eigenvector with an eigenvalue e. Then the eigenvalue 
problem reduces to the following difference equation 



ea r , 



t±(a n -i + a n+ i), 



(2-3) 



with the boundary condition ao = «at+i = 0. Assuming a n ~ e , it can be shown 
thatED 



2 1± cos 6* r , 



a, 



N+l 



(sin r , sin 26 r , ■ ■ ■ , sin M6 r ), 



(2-4) 



where r = l,2,...,iVis the chain eigenvalue index and 9 r = rir/(N + 1). Note that 
odd N chains have a zero-energy eigenstate with an eigenvector that has nonzero 
amplitudes, constant in magnitude and alternating in sign, on the sublattice of the 
chain ends. 
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2.2.2. AA stacking 

Although AA stacking is not energetically favorable, it is still interesting to 
consider this arrangement for pedagogical purposes. In the case of AA stacking, the 
Hamiltonian at K is given by 



#aa(p) 
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(2-5) 
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where tt = p x + ip y . 

As we now explain, the electronic structure of AA stacked iV-layer graphene 
can be thought of as consisting of separate ID chains for each wavevector in the 
2D triangular lattice Brillouin zone of a single graphene layer. For an eigenvector 
(ai,b\,--- ,aAr,67v) with an eigenvalue e and fixed 2D momentum, the difference 
equations in this case are 



ea n = t±(a n -i + a n+ i) + U7r'6 n , 
eb n = t±(b n -x + K+i) + vTTa n , 



(2-6) 



with the boundary condition ao = aAr + i = bo = bjy+i = 0. 

Let c n = a n + b n e~ 1 ^ and d n = a n — b n e~ l< ^ where cj) = tan -1 (p y /p x ), then 



(e - v\p\)Cn = t±(c n -l + Cn+l), 

(e + v\p\)d n = t±(d n -i + dn+i), 



(2-7) 



with the same boundary condition co = cn+i = do = d]y + i = 0. Thus the energy 
spectrum is given by 



±u|p| + 2tj_ cos 



rir 



N + 1 



(2-f 



where r = 1, 2, • • • , N. Note that for odd N, the r = (iV + l)/2 mode provides two 
zero-energy states at p = 0. 

Figure [3] shows the band structure of AA stacked trilayer and tetralayer graphene 
near the K point. Because of the hybridization between a-a and (3-(5 sublattices in 
each layer, zero-energy states occur at momenta that are remote from the K and K' 
points. In the following we turn our attention to stacks in which adjacent graphene 
layers have a relative rotation of 60 degrees. As we show, in this case the zero-energy 
states always occur precisely at the Brillouin-zone corners. 
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Fig. 3. Band structure near the K point for (a) trilayer and (b) tetralayer graphene with AA stack- 
ing for nearest intralayer neighbor hopping t = 3 eV and nearest interlayer neighbor hopping 
t± = O.lt. 



2.2.3. AB stacking 

In the case of AB stacking, the Hamiltonian at K has the following form, 



H AB (p) 



( o 

vtt 
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(2-9) 



We will see that the subtle difference in the Hamiltonian compared to the AA case 
changes the electronic structure in a qualitative way. To obtain the energy spec- 
trum of AB stacked iV-layer graphene, let us consider corresponding difference equa- 



tions 



ID 



ea 2n -i = (W)6 2 n-i, 

£&2ri-l = i±(a2n-2 + 02n) + (wr)a2n-l, 

ea 2n = t ± (b 2n -i + fon+l) + {vTT ] )b 2n , 

eb 2n = {vTr)a 2n , (2-10) 

with the boundary condition ao = ajy+i = bo = bjy+i = 0. 

Let c 2n -\ = b 2n -i and c 2n = a 2n , then the difference equations reduce to 

(e - v 2 \p\ 2 /e)c n = t±(Cn-l + Cn+l), (2-11) 

with the boundary condition cq = cm+i = 0. Then the energy spectrum is given by 



e — v 2 \p\ 2 je = 2t± cos 



TTT 

N + 1 



(2-12) 
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Fig. 4. Band structure near the K point for (a) trilayer and (b) tetralayer graphene with AB stack- 
ing for nearest intralayer neighbor hopping t = 3 eV and nearest interlayer neighbor hopping 
t± = O.lt. 



where r = 1, 2, • • • ,N. Thus 



< P = ^ - [wtij ± r v + tl cos2 [wh ) ■ (2 - i3) 

Note that relativistic energy spectrum for a particle with the momentum p and 
mass m is given by 

e p = ^\p\ 2 c 2 + m 2 c i . (2-14) 



Thus we can identify m r v 2 = t± cos 

For a massive mode with mass m r , the 



as the effective mass for mode r. 
ow-energy spectrum is given by 




(2-15) 



For odd N, the mode with r = (JV + l)/2 is massless and its energy is given by 

£±&±v\p\. (2-16) 

For even N, all N modes are massive at low energies. Therefore, the low-energy 
spectrum with odd number of layers is a combination of one massless Dirac mode 
and N—l massive Dirac modes, whereas the low-energy spectrum with even number 
of layers is composed of only massive Dirac modes. 

Figure [4] shows the band structure of AB stacked trilayer and tetralayer graphene 
near the K point. As discussed earlier, the trilayer has one massless mode and two 
massive modes, while the tetralayer has all massive modes at low energies. Note that 
at p = 0, each massless mode gives two zero energies while each massive mode gives 
one zero energy. Therefore, for odd N, there are 2 + (N — 1) = N + 1 zero-energy 
states while for even N, there are iV zero-energy states. 



Electronic Structure of Multilayer Graphene 



7 



2.2.4. ABC stacking 

In the case of ABC stacking, the Hamiltonian at K is given by 

\ 
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(2-17) 
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Unfortunately for ABC stacking, there do not exist low-order difference equations 
with a simple boundary condition. Instead we can easily derive a low-energy effec- 
tive Hamiltonian. Surprisingly, it turns out that ABC stacked TV-layer graphene is 
described by iV-chiral 2D electron system. (More detailed discussion for the effective 
theory of arbitrarily stacked graphene will be presented in fj3]) 

It is important to recognize that in ABC stacking, there is vertical hopping 
between all the lower layer (5 sites and all the upper layer a sites. For ir = each 
a- (5 pair forms a symmetric-antisymmetric doublet with energies ±ij_, leaving the 
bottom ot\ and top (5^ sites as the only low-energy states. This behavior is readily 
understood from the stacking diagrams, in Fig.[2j It is possible to construct a 2 x 2 ir- 
dependent low-energy effective Hamiltonian for the low-energy part of the spectrum 
using perturbation theory. The same procedure can then be extended to arbitrary 
stacking sequences. 

The simplest example is bilayer graphene, 15 ' Low and high energy subspaces 
are identified by finding the spectrum at ir = and identifying all the zero-energy 
eigenstates. The intralayer tunneling term, which is proportional to it, couples low 
and high energy states. Using degenerate state perturbation theory, the effective 
Hamiltonian in the low energy space is given to leading (2nd) order in it by 



H e 2 ff iP) 




Ml! 

2m 




-t± 







(„t)* 





(2-18) 



where we have used a (ai,/?2) basis, m = t±/2v 2 and v = vir/t±. In the same way 
we find that the effective Hamiltonian of ABC stacked iV-layer graphene is given by 



H e J f (p) 



-t± 







(„t)" 




(2-19) 



using a (ai,/3jv) basis. The leading correction appears at order N in it because the 
unperturbed high-energy states are localized on a ai+i) pair and the perturbation 
is intralayer tunneling. Note that we have for mathematical convenience chosen a 
gauge in which the single-layer Hamiltonian is 



(P) 



VTT 
VTT 



(2-20) 
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We can prove Eq. (2-19) by the mathematical induction method. Imagine that 
we add one more layer on top of iV-layer graphene with ABC stacking. Then the 
combined Hamiltonian is given by 



H eff 



1(P) 



-t± 



( 





(„t)tf 
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(2-21) 



using a (a±, /3jv> ajv+i, Pn+i) basis. 

Let P be a low-energy subspace spanned by (ai,/3jv+i) and Q be a high-energy 
subspace spanned by (ajv+ij Pn)- Note that the effective Hamiltonian can be derived 
using the degenerate state perturbation theory^ 



H, 



eff ~ Hpp — HpQ 



1 



if, 



-H, 



QQ 



QP- 



(2-22) 



Here the projected Hamiltonian matrices to P and Q subspace are given by 

H QQ {p) = ti ( J J ) , ^PQ(P) = -*± ( ° ^ ) , (2-23) 
and Hpp{p) = 0. Thus we can easily show that, 

" L (P) - -'i 







(l/t)"+l 





(2-24) 



which proves Eq. (2-19). The corresponding energy spectrum in Eq. (2-19) is given 

by 

' V \P\ 



'eff,P 



±t± 



t± 



N 



(2-25) 



Figure [5] shows the band structure of ABC stacked trilayer and tetralayer graphene 
near the K point. Note that at p = 0, there are only two zero energy states no matter 
how thick the stack is. 

2.2.5. Arbitrary stacking 

It is easy to generalize the previous discussion to construct the Hamiltonian for 
an arbitrarily stacked multilayer graphene system. Figure [6] shows the band structure 
of ABCB stacked tetralayer graphene and ABBC stacked tetralayer graphene near 
the K point. For the ABCB stacked tetralayer graphene, the low-energy spectrum 
looks like a superposition of a linear dispersion and a cubic one. For the ABBA 
stacked tetralayer graphene, zero energies appear not only at the Dirac point but 
also away from it. A more detailed low-energy spectrum analysis will be presented 
in© 

2.3. Landau level spectrum 
2.3.1. Preliminaries 

In the presence of a magnetic field B 

P ~ 



p + ~A, where A is the vector potential with B 



Bz, a Hamiltonian is modified by 
V x A. The quantum 
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Fig. 5. Band structure near the K point for (a) trilayer and (b) tetralayer graphene with ABC 
stacking for nearest intralayer neighbor hopping t — 3 eV and nearest interlayer neighbor hop- 
ping t± = O.li. 




Fig. 6. Band structure near the K point for tetralayer graphene with (a) ABCB stacking and (b) 
ABBC stacking for nearest intralayer neighbor hopping t = 3 eV and nearest interlayer neighbor 
hopping t± — O.lt. 



Hamiltonian is most easily diagonalized by introducing raising and lowering opera- 
tors, a = lifi I\f2% and o) = ln/\/2h substitution, where I = y/hc/e\B\, and noting 
that [a, at] = 1. We can then expand the wavefunction amplitude on each sublattice 
of each layer in terms of parabolic band Landau level states | n) which are eigenstates 
of the a)a. For many Hamiltonians, including those studied here, the Hamiltonian 
can be block diagonalized by fixing the parabolic band Landau-level offset between 
different sublattices and between different layers. This procedure is familiar from 
theories of Landau-level structure in other multiband k ■ p theories. 
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Fig. 7. Landau levels of (a) trilayer and (b) tetralayer graphene with AA stacking for nearest 
intralayer neighbor hopping t = 3 eV and nearest interlayer neighbor hopping t± — O.lt. Landau 
levels were shown up to n = 10. 



2.3.2. A A stacking 

In the case of AA stacking, let us choose the n-th Landau level basis at K as 
(ai, n -i,/3i,n, • • • ,aN,n-i,PN,n)- Then Eq. reduces to 



HAA{n) 
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(2-26) 
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where e n = \/2nhv/l. Note that 2D Landau level states with a negative index do not 
exist so the corresponding basis states and matrix elements are understood as being 
absent in the matrix block. Thus H^^n = 0) is a N x N matrix, while Hj^\(n > 0) 
is a 2iV x 2N matrix. 



By diagonalizing Eq. (2-26) using the difference equation method, we can obtain 



the exact Landau level spectrum. For n > 0, Landau levels are given by 



±e n + 2ij_ cos 



T7T 



N+l 



(2-27) 



where r 
2t± cos (- 



1,2, ••■ ,N. Note that for n = 0, Landau levels are given by e T) q = 
N _^lj- Thus for odd N, there exists one (^-independent) zero-energy Lan- 
dau level at r = (N + l)/2. 

Figure [7] shows the Landau levels of AA stacked trilayer and tetralayer graphene 
as a function of magnetic fields. For the trilayer, there is one zero-energy Landau 
level, while for the tetralayer, there is no zero-energy Landau level. Note that there 
are Landau levels crossing the zero-energy line in AA stacking. 
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2.3.3. AB stacking 

In the case of AB stacking, a proper choice of the n-th Landau level basis at K 
is (ai, n _i, Pi >n , ot2,n, Ih,n+l,ac3,n-l,{k,n, «4,n, Pi,n+i, ■••) such that all the interlayer 



hopping terms are contained in the n-th Landau level Hamiltonian. Then Eq. (2-9) 
reduces to 
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(2-28) 



/ 



where e n = y/2nhv/l. As discussed earlier, special care should be given for states 
with a negative index. 



For the Hamiltonian in Eq. (2-28), there do not exist corresponding difference 



equations with a proper boundary condition, thus cannot be diagonalized analyti- 



cally. From Eq. (2-15 ), however, we can find the low-energy Landau levels for massive 
mode with mass m r as 



+hoj r ^Jn{n + 1) if t± cos 



-huj r \/n{n + 1) if t± cos ( -jP^ ) > 0, 




<0, 



(2-29) 



where uj r = e\B\/m r c and r = 1,2, ■•• ,N, which is proportional to B. These 
equations apply at small B, just as the B = limiting low-energy dispersions applied 



at small momentum ir. For the massless mode, from Eq. (2-16) Landau levels are 
given by 

et = ±e„, (2-30) 

which is proportional to B 1 / 2 . 

Figure [8] shows the Landau levels of AB stacked trilayer and tetralayer graphene 
as a function of magnetic fields. Note that the linear B dependence expected for 
massive modes applies over a more limited field range when the mass is small. For 
the trilayer, Landau levels are composed of massless Dirac spectra (oc i? 1 / 2 ) and 
massive Dirac spectra (oc B), while for the tetralayer, Landau levels are all massive 
Dirac spectra. This is consistent with the band structure analysis shown in Fig. |4j 

Note that the massive modes in Eq. (2-29) have two zero-energy Landau levels 
for n = — 1 and 0, whereas the massless mode in Eq. (2-30) has one for n = 0. There 
are therefore N zero-energy Landau levels in both even and odd iV AB stacks. This 
property can also be understood directly from the Hamiltonian in Eq. (2-28), by 
eliminating negative n basis states and rearranging rows to block diagonalize the 
matrix. 

2.3.4. ABC stacking 

In the case of ABC stacking, a proper choice of the n-th Landau level basis at K is 
(ai >n _i, #L,n, at2, n , (h,n+i, a3,n+i > Ps,n+2, ■ ■ ■ ) such that all the interlayer hopping terms 
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are contained in the n-th Landau level Hamiltonian. Then Eq. (2-17) reduces to 

\ 
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(2-31) 
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where e n = yj2nhv/l. 

The low-energy spectrum can be obtained from the effective Hamiltonian in 



Eq. (2-17). For n > 0, Landau levels are given by 

e± = ±HuJ NV / n(n + 1) • • • (n + N - 1) 



(2-32) 



where = t±(V2hv/t±l) N oc fi^/ 2 , while for n = —N + 1, —N + 2, • • • ,0 they 
are zero. Note that there are iV zero-energy Landau levels for ABC stacked iV-layer 
graphene. 

Figure|9]shows the Landau levels of ABC stacked trilayer and tetralayer graphene 
as a function of magnetic fields. For the trilayer, Landau levels are proportional to 
.B 3 / 2 , while for the tetralayer, Landau levels are proportional to B 2 . 

2.3.5. Arbitrary stacking 

It is straightforward to generalize the previous discussion to construct the Hamil- 
tonian in Landau level basis for an arbitrarily stacked multilayer graphene system. 
As discussed earlier, special care should be given for states with a negative index. 
Figure [To] shows Landau levels of ABCB stacked tetralayer graphene and ABBC 
stacked tetralayer graphene. For the ABCB stacked tetralayer graphene, the Lan- 
dau levels look like a superposition of B 1 / 2 and B 3 / 2 levels, which is consistent with 
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Fig.[6^a). For the ABBA stacked tetralayer graphene, there are Landau levels cross- 
ing the zero-energy line, which is consistent with Fig. |6^b). Detailed low-energy 
Landau-level spectrum analysis will be presented in ^3} 

2.4. Quantum Hall conductivity 

Applying the Kubo formula to a disorder-free systems we find that the conduc- 
tivity tensor with an external magnetic field along z is given by 



(2-33) 
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where f n is Fermi factor of n-th energy state, i,j = x, y and 

(n\ hvi \m) (m\ hvj |n) (m\ hvi \n) (n\ hvj \m) 



{e n - e m )(e n -e m + huj + irf) (e n - e m )(e n - e m - Huj - irj) 

(2-34) 

Here Vi is a velocity operator obtained by taking a derivative of the Hamiltonian 
H{p) with respect to pi. Note that in case of multilayer graphene, the velocity 
operator is constant, i.e. it does not depend on the Landau level index. 

The appropriate quantized Hall conductivity is obtained by evaluating an = 
a xy (0). Detailed analysis of the quantum Hall conductivity will be presented in 

§3. Chiral decomposition of energy spectrum 



In this sectiorp] we demonstrate an unanticipated low-energy property of 
graphene multilayers, which follows from an interplay between interlayer tunneling 
and the chiral properties of low-energy quasiparticles in an isolated graphene sheet. 
Our conclusions apply in the strongest form to models with only nearest-neighbor 
interlayer tunneling, but are valid over a broad field range as we explain below. We 
find that the low-energy band structure of any graphene multilayer consists of a set 
of independent pseudospin doublets. Within each doublet, the bands are described 
by a pseudospin Hamiltonian of the form 

Hj(k) oc k J [cos{J0 k )T x ± sm(J0 k )T y ], (3-1) 

where T a is a Pauli matrix acting on the doublet pseudospin, k is an envelope 
function momentum measured from either the K or K' corner of the honeycomb 
lattice's Brillouin-zone|3'[ll' k = \k\, and is the orientation of k. The ± sign in 



Eq. (3-1) assumes the opposite signs in graphene's K and K' valleys. Following the 
earlier work on graphene bilayersji^l we refer to J as the chirality index of a doublet. 
In the presence of a perpendicular magnetic field B, Hj{k) yields J Landau levels 
at E = and i? / levels with \E\ oc B J I 2 . Taking the twofold spin and valley 
degeneracies into account, the number of independent zero-energy band eigenstates 
at the Dirac point (fc = 0) is therefore 8iV£), where Nd is the number of pseudospin 
doublets. We find that, although Nd depends on the details of the stacking sequence, 

N D 

Y^Ji=N (3-2) 

i=l 



in an N- layer stack. It follows from Eq. (3-2 ) that the Hall conductivity of an iV-layer 



stack has strong integer quantum Hall effects with plateau conductivities, 

4e 2 (N \ 
<7xy = ±— ( — + n) , (3-3) 



where n is a non-negative integer. 



*' The content of this section provides a more complete explanation of the arguments presented 
earlier in Ref. 10). 
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3.1. Partitioning rules 

The low-energy band and the Landau level structure can be read off the stacking 
diagrams illustrated in Fig.[2]by partitioning a stack using the following rules, which 
are justified in the following section. 

(i) Identify the longest nonoverlapping segments within which there are no re- 
versals of stacking sense. When there is ambiguity in the selection of nonoverlapping 
segments, choose the partitioning which incorporates the largest number of layers. 
Each segment (including for interior segments the end layers at which reversals take 
place) defines a J-layer partition of the stack and may be associated with a chirality 
J doublet. 

(ii) Iteratively partition the remaining segments of the stack into smaller J 
elements, excluding layers contained within previously identified partitions, until all 
layers are exhausted. 

The chirality decompositions which follow from these rules are summarized in 
Table |IJ Note that this procedure can result in J = 1 doublets associated with 
separated single layers which remain at the last step in the partitioning process. 

In applying these rules, the simplest case is cyclic ABC stacking for which there 
are no stacking sense reversals and therefore a single J = N partition. In the 
opposite limit, AB stacking, the stacking sense is reversed in every layer and the 
rules imply N/2 partitions with J = 2 for even N, and when N is odd a remaining 
J = 1 partition. Between these two limits, a rich variety of qualitatively distinct 
low-energy behaviors occur. For example, in the ABCB stacked tetralayer, ABC is 
identified as a J = 3 doublet and the remaining B layer gives a J = 1 doublet. The 
low-energy band structure and the Landau level structure of this stack, as illustrated 
in Figs. [6^a) and 10 [a), have two sets of low-energy bands with \E\ oc k, k 3 , Landau 



levels with \E\ oc iF/ 2 , £? 3 / 2 , and four zero-energy Landau levels per spin and valley. 
All these properties are predicted by the partitioning rules. We have explicitly 
checked that the rules correctly reproduce the low-energy electronic structure for all 
stacking sequences up to N = 7. Because each layer is a member of one and only 



one partition, the partitioning rules imply the chirality sum rule in Eq. (3-2). 

3.2. Degenerate state perturbation theory 

We start from the well-known J = 1 massless Dirac equatiorI3'II) k p model for 
isolated sheets, 

"o)< < 34 > 

where ir = p x + ip y and v is the quasiparticle velocity. In the presence of an external 
magnetic field, it and 7T are proportional to the Landau level raising and lowering 



operators, so that Eq. (3-4) implies the presence of one macroscopically degenerate 
Landau level at the Dirac point for e ach spin and valley, and therefore, to the N = 1 
quantum Hall effect®® of Eq. (p-3b. An iV-layer stack has a two-dimensional band 



structure with 2N atoms per unit cell. The Hamiltonian can be written as 



H — H± + , 



(3-5) 
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Table I. Chirality decomposition for N = 3, 4, 5, 6 layer stacks. 
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ABCBC 


3+2 


ABABAB 


2+2+2 


ABCBA 
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1+4+1 


ABACB 


1+4 


ABACBA 


1+5 



where H± accounts for interlayer tunneling and Hu for intralayer tunneling. Hii 
is the direct product of massless Dirac model Hamiltonians Hmd for the sublattice 
pseudospin degrees of freedom of each layer. We construct a low-energy Hamiltonian 
by first identifying the zero-energy eigenstates of H± and then treating Hu as a 
perturbation. 

Referring to Fig. [2] we see that H± is the direct product of a set of finite-length 
ID tight-binding chains, as shown in Eq. (2-4), and a null matrix with dimension 
equal to the number of isolated sites. The set of zero-energy eigenstates of H± 
consists of the states localized on isolated sites and the single zero-energy eigenstates 
of each odd-length chain. 

The low-energy effective Hamiltonian is evaluated by applying leading order de- 
generate state perturbation theory to the zero-energy subspace. The matrix element 
of_the effective Hamiltonian between degenerate zero-energy states r and r' is given 
by 



{&r\H\V r ,) = {V r \H\\ Q{-Hl l )QH\ 



n—l 



(3-6) 



where n is the smallest positive integer for which the matrix element is nonzero, and 
Q = 1 — P, P is a projection operator onto the zero-energy subspace. To understand 
the structure of this Hamiltonian, it is helpful to start with some simple examples. 

3.2.1. ABC stacking 

For ABC stacked AMayer graphene, the zero-energy states are the two isolated 
site states in bottom and top layers, a\ and /3n- N — 1 sets of two-site chains form 
high-energy states. Because H\\ is diagonal in layer index and H± (and hence Hj}) 
can change the layer index by one unit, the lowest order at which a± and /3jy are 
coupled is n = N. 
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According to Eq. (2-4), the wavefunction of each two-site chain is given by 



1 



V2 



( \(3 r ) + a r \a r+ i)), 



(3-7) 



with the energy e r = t±a r , where a r = ±1 and r = 1, 2, • • • , N — 1. From Eq. (3-6 ), 



{<*i\H\p N ) 



N-l 



(ai| H { \ [Qi-H^QH^ \p N ) 



(-£l) • • • (-EjV-l) 



, (-o"i/2) • • • (-ctjv-i/2) . t jv 
"'^ (-ai)-(-^-i) ^ 



{oy} 



UN 



1 



)AT-1 



01,— jCTjv-I 



-^(- t ) 7V , 



(3-8) 



where i/ = vir/t±. Here (ai|V|<2> CTl ) = -(l/V2)t ± i/t, <^ CTJV _ 1 1 ^|/?at> = -((Tjv-i 
/V^)t±^ and (^ayl^l^ay+i) = -(cr r /2)4_L^ were used. Thus, the effective Hamil- 
tonian of iV-layer graphene with ABC stacking has a single J = N doublet given 
by 

(i/t)* 



e// 
TV 



iV 







(3-9) 



3.2.2. AB stacking 

For AB stacked iV-layer graphene, the high-energy Hilbert space consists of a 
single iV-site ID chain, excluding its zero-energy eigenstate when N is odd. There 
is an isolated site in each layer which is connected to both its neighbors at or- 
der n = 2 forming an isolated site chain. When N is even, this chain is diag- 
onalized by N/2, J = 2 doublets formed between a-sublattice and /5-sublattice 
chain states0'O3'[^lJ'[iH]' When N is odd, the zero-energy chain state is mapped 
to an equal-magnitude oscillating-sign linear combination of isolated site states by 
intralayer tunneling at order n = 1, yielding a J = 1 doublet. The {N — l)/2, J = 2 
doublets are then formed between a-sublattice and /3-sublattice isolated site chain 
states in the orthogonal portion of the isolated state subspace. 

Let us consider the low-energy spectrum of AB stacking in more detail. From 



Eq. (2-4) energy spectra and wavefunctions of the single iV-site chain are given by 



e r = 2t± cos 9 r , 



12V 



N + 



-( sin0 r + sin 26 r \a 2 ) + sin 2,0 r \/3 3 ) + sm40 r |a 4 ) • • • ) , (3-10) 



where 6 r = -^-j- and r = 1, 2, • • • ,N. 

First, let us consider the case with even N. Then the low-energy states come 
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from the isolated sites or equivalently their superpositions. Let us define 



\9 r ) = y_^_( s in^ e -^ |ai)+sin20 r e^ |£ 2 )+sin30 r e-^ |a 3 )+sin40 r e^ |/3 4 ) • • •) 

(3-11) 

such that 

(&r\V\$ s ) = -S r ,sW\, (3-12) 

where e*^ = Then the matrix elements between the low-energy states are 

given by the second order perturbation theory: 

mHm _ £ Mgj _ _ V((i/£r)M2 . (3.13) 

s—1 

Note that e r = —en+i-t and these two modes form a 2-chiral system with energies 
±|e r |. The chirality can be manifested clearly if we define 

e i<t> 

l«r) = ^ (l«V) + 1%+1-r}) , 

A-) = ^-(|!?r)-|S^+l-r». (3-14) 
Then the Hamiltonian of the 2-chiral system for r = 1, 2, • • • , N/2 is given by 



in a (a r ,/3 r ) basis with m r u = tj_ cos ^^qr[J. Thus the system is described by a 

combination of massive Dirac modes with different masses. 

For odd N, there is a zero-energy state in the iV-site chain at r = (JV + l)/2 in 
Eq. (3-10). Thus in addition to the massive modes, there exists one massless Dirac 
mode, 

/ <Pn+i\V\<Pn+i) = -\v\. (3-16) 



2 2 

Thus the system is described by one massless Dirac mode and a combination of 
massive Dirac modes with different masses. 

3.2.3. ABC+B type stacking 

A more complex and more typical example is realized by placing a single reversed 
layer on top of ABC stacked iV-layer graphene with N > 2. Note that the last chain 
has three sites, thus it has a zero-energy state defined by 

\Pn+i) = ^(\Pn+i)-\Pn-i)), (3-17) 
and two high-energy states defined by 

I*™-!) = ^= (Wi> + ^ \«n) + , (3-18) 
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where cr^-i = ±1- Then the first-order perturbation theory gives 

t± 



otN+i\H\j3 r 



F"t, 



N+l/ ~ y/2 

suggesting the existence of the massless Dirac mode with reduced velocity. 



(3-19) 



Similarly as Eq. (3-8), we obtain 



H eff 



V 











(„t)2 
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sft 
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V2 



















V2 


2 




» 2 

2 





2 





/ 



(3-20) 



using a (ajv+i, Pn+i: a i> Pn) basis. The first 2x2 block in Eq. (3-20) gives a J = 1 
doublet with a reduced velocity. Note that the matrix in Eq. (3-20) is not block 
diagonal thus we cannot simply say that the second 2x2 matrix block is a N- 
chiral system. The J = N doublet in this instance includes both the (ai,/3j\r) 
subspace contribution and an equal contribution due to perturbative coupling to th e 
(on+i, (3^f +1 ) subspace. Using a similar perturbation theory shown in Eq. (2-22), 



we can obtain higher order correction by integrating out the massless Dirac mode 
which forms a higher energy state. Then the final Hamiltonian is reduced to 



H eff 



Hi ® H N , 



where 



-t± 





v/V2 







H 



N 



-t± 







(„t)" 




(3-21) 



(3-22) 



This means that the combined system can be described by a combination of one 
1-chiral system with reduced velocity and one iV-chiral system. Note that stacking 
a layer with an opposite handedness partitions a system into systems with different 
chiralities. 

3.2.4. Arbitrary stacking 

The relationship between the electronic structure of a general stack and the 
partitioning procedure explained above can be understood as follows. 

(i) First, note that a partition with chirality J has isolated sites in its terminal 
layers that are coupled at order J in perturbation theory. In the case of J = 1 
partition, the chain opposite to the single isolated site always has an odd length and 
provides the zero-energy partner; isolated site to chain coupling therefore always 
occurs at first order. 

(ii) Next, consider the perturbation theory, truncating at successively higher 
orders. When truncated at first order, the J = 1 partitions are isolated by higher 
J blocks within which the Hamiltonian vanishes. Each J = 1 partition therefore 
yields a separate massless Dirac equation with velocitiea*^! that can be smaller than 



*' The velocity of the J = 1 doublets is determined by the strength of the coupling between 
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the graphene sheet Dirac velocity. When the perturbation theory is truncated at 
second order, the Hamiltonian becomes nonzero within the J = 2 partitions. The 
eigenenergies within the J = 1 partitions are parametrically larger, and the Hamil- 
tonian within the J > 2 partitions is still zero. To leading order therefore, the J = 2 
partitions are separated, and their isolated states are coupled at the second order in 
perturbation theory so that each provides a J = 2 doublet such as that of an isolated 
bilayer. If two or more J = 2 partitions are adjacent, then their Hamiltonians do not 
separate. In this case, there is a chain of second order couplings between isolated 
states, such as those of an even-length AB stack, but the end result is still J = 2 
doublet for each J = 2 partition. 

(iii) The identification between partitions and chiral doublets can be established 
by continuing this consideration up to the highest values of J which occur for a 
particular stack. 

(iv) Then, the effective Hamiltonian of any iV-layer graphene is as follows: 

H e J f « * ® fl> a ® • • • ® Hj Nd , (3-23) 



with the sum rule in Eq. (3-2). Note that Np is half the sum of the number of 
isolated sites and the number of odd-length chains. 

3.3. Discussion 

3.3.1. Effects of remote hopping 

The minimal model we have used to derive these results is approximately valid 
in the broad intermediate magnetic field B range between ~ 10 and ~ 100 T, 
over which the intralayer hopping energy in field (~ hv/i where i = \Jhcfe\B\ ~ 
25 nm/[B(T)] 1 / 2 is the magnetic length) is larger than the distant neighbor interlayer 
hopping amplitudes that we have neglected (72 ~ —20 meV), but still smaller than 

For example, if we consider a\ — > hopping process in ABA stacked trilayer in 
Fig. [2] the valid range of magnetic field for the minimal model is given by 

(hv/l) 2 

I72I < < t±- (3-24) 

t± 

When 72 does not play an important role (in N = 2 stacks, for example), the 
lower limit of the validity range is parametrically smaller. The minimum field in 
bilayers has been estimated to be ~ 1 TJI5I by comparing intralayer hopping with 
the 73 ~ 0.3 eV interlayer hopping amplitude, 

hv 3 /l < ^-/-L < t ± , (3-25) 

where 1*3 = (V3/2)ajs/h and a is a lattice constant of graphene. 



the odd-length chain zero-energy state and isolated states on the sublattice opposite to the chain 
ends. For a chain of length 2N — 1, the chain's zero-energy state has nonzero amplitude on the TV 
odd-index sites. The velocity is reduced from the single sheet velocity by a factor of v/ M/N, where 
M is the number of isolated sites opposite to the N odd-index sites. In a similar manner, higher J 
doublet Hamiltonians are sometimes altered by a multiplicative factor by perturbative coupling to 
smaller J doublets as in the single reversed layer example. 
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Fig. 12. Landau level spectrum near the K valley as a function of 72 for an ABA stacked trilayer for 
(a) B = 1 T and (b) B = 10 T. Here t = 3 eV, t± = O.lt, and w c = eB/mc, with m = t±/2v 2 , 
were used. Note that for this case the Landau level structures near K and K' valleys are not 
identical. 



Figures 11 and |12| show the Landau level spectrum at the K valley as a function 
of 73 for an AB stacked bilayer, and as a function of 72 for an ABA stacked trilayer, 
respectively. In the case of the bilayer, the dependence of the Landau levels on 73 
is weak for B larger than 1 T, whereas in the case of the trilayer, the Landau level 
spectrum still strongly depends on 72 for B = 1 T, but the dependence becomes 
weak for B above 10 T, confirming the above argument. 
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3.3.2. Quantum Hall effect 
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Fig. 13. (Color online) Noninteracting system Hall conductivity as a function of the Fermi energy 
for all inequivalent four-layer graphene stacks when B = 10 T, t = 3 eV, and t± = O.lt. The 
dependence of the Hall conductivity on Fermi energy is simply related to the dependence of the 
Hall conductivity on total electron density. The Hall conductivity calculations shown in this 
figure assume neutralizing ionized donors spread equally between the four layers. 



In Fig. 13 we plot the noninteracting Hall conductivity as a function of Fermi 
energy for the four distinct four-layer stacks. When electron-electron interactions 
are included at an electrostatic mean-field (Hartree) level and the neutralizing ion- 
ized dopants (responsible for the Fermi energy shift away from the Dirac point) are 
assumed to be equally distributed among the layers, the Landau levels with E ^ 
are shifted by electrostatic potential differences between the layers. There is, how- 
ever, no influence of electrostatics on the E = levels. This property follows from 
the perfect particle-hole symmetry of the models we employ, which implies a uni- 
form charge distribution among the layers at the neutrality point. Remote (72 2nd 
neighbor) interlayer hopping does shift the E = Landau level in the ABAB stacked 
tetralayer and weakly lifts the degeneracy responsible for the large jump between 
the ±(4e 2 /h)N/2 Hall plateaus. This example demonstrates a tendency toward the 
grouping of N spin and valley degenerate Landau levels very close to E = in 
general iV-layer stacks even when remote neighbor hopping is included. Small gaps 
between these Landau levels are unlikely to lead to Hall plateaus unless disorder is 
very weak. When disorder is weak, on the other hand, electron-electron interaction 
effects beyond Hartree level are likely to be important and lead to strong quantum 
Hall effects at many filling factors, often ones associated with broken symmetries of 
different types 113 @^E§>IIi}'[Mf The property that the Hall conductivity will tend 
to jump by four units on crossing the Dirac point for arbitrarily stacked tetralayer 



Electronic Structure of Multilayer Graphene 



23 



(a) a 4 ©0 4 (b) a 4 @j3 4 





a. 


>® 




• 




f ® 


03 






,® 




> 




>® 


B 2 



Fig. 14. (Color online) Stacking diagrams for tetralayer graphene with (a) ABBC stacking and (b) 
ABBA stacking. Shaded ovals link nearest interlayer neighbors. 



graphene is the most obvious experimental manifestation of the chirality sum rule 
discussed in this paper. In practice charged multilayers (Ep ^ 0) would normally 
be prepared by placing the system on one side of an electrode and gating. Even 
though gating will redistribute charge and shift electric potentials differently in dif- 
ferent layers, the Landau level bunching we discussed should still be clearly reflected 
in quantum Hall effect measurements. 

3.3.3. Effects of the same stacking inside 

The analysis presented so far is based on the assumption that stacking one layer 
directly on top of its neighbor, AA stacking, is not allowed. When interior AA 
stacking does occur, we can still apply a similar diagram analysis and identify the 
zero-energy states at the Dirac point. In this case, however, zero-energy states can 
appear not only at the Dirac points but also at other points in momentum space. The 
degenerate state perturbation theory at the Dirac point discussed so far therefore 
cannot completely capture the low-energy states. 

As an example, let us consider ABBC stacked tetralayer graphene, as illustrated 
in Fig. |14[a). Here, in addition to oc\ and f3 4 , there are two zero-energy states at 
each three-site-chain defined by 



&) = 4O01>-|O3» 



V2 

|« 4 ) = _L(| a4 )_| /?2 )). ( 3 .26) 

Thus the matrix elements between low-energy states are given by 

ai|ff|&) = {a 4 \H\(3 4 ) = ~^jf ] - (3-27) 

Therefore the system can be described by two massless Dirac modes with reduced 
velocity, as shown in Figs.^b) and 10 [b). 



Another example is ABBA stacked tetralayer graphene, as illustrated in Fig. 14 (b) 



In this case, there are two zero-energy states at a± and a 4 . The high-energy states 



<P r are given by Eq. (2-4) with N = A, thus we get 

4 



/ I HI \ V {ai\V\$r) {$r\V\(X4) , , ,2 

(a:i|-n \a 4 ) = y j r = —ct\_\v\ , (3-28) 

1 \ e r) 
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where c = g J2 r sin (~) sin (^p) /cos (t?) = —1. Here the low-energy state is 
composed of one non-chiral massive mode. Note that because of the non-chirality, 
there are no zero-energy Landau levels. 

3.3.4. Pseudospin magnetism 
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Fig. 15. (Color online) In-plane projected pseudospin orientation of (a) J = 1, (b) J = 2, (c) J = 3 
and (d) J = 4 chiral 2D electron system for a neutral, unbiased system with coupling constant 
a = e 2 /ehv = 1 where e is the dielectric constant. For J > 1, the arrows are shorter in the core 
of the momentum space vortex because the pseudospins in the core have rotated spontaneously 
toward z or — z direction indicating the pseudospin magnetic state. 



Finally, we note that in the presence of electron-electron interactions, chiral two- 
dimensional electron system (C2DES) tends toward momentum-space vortex states 
in which charge is spontaneously shifted between layers^ and that these instabilities 
are stronger in systems with larger J. Figure 15 shows in-plane projected pseudospin 
orientation for J = 1, 2, 3, 4 C2DESs, which correspond to N = 1, 2, 3, 4 ABC-stacked 
graphene multilayers. Note that for J > 1, the arrows in the core of the momentum 
space vortex have rotated spontaneously toward z or — z direction indicating the 
spontaneous charge transfer between layers. The present work identifies ABC stacked 
multilayer graphene as the most likely candidate for this particular type of exotic 
broken symmetry state. Other types of broken symmetry might occur for other 
stacking sequences, especially in the quantum Hall regime. 
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§4. Conclusions 



We have shown that iV-layer graphene at intermediate and strong magnetic fields 
has a strong tendency towards the appearance of N spin and orbitally degenerate 
Landau levels very close to E = 0. This property should lead to strong quantum 
Hall effects at ±(4e 2 /h)N/2 in many iV-layer stacks. The origin of this behavior is 
the following chirality sum rule: i) The low-energy bands of multilayer graphene can 
be decomposed into Nd doublets with chirality Jj. ii) Although Np depends on the 
stacking sequence, X^i=a Ji = N in an iV-layer stack. 

The chirality sum rule applies precisely only to idealized models with only 
nearest-neighbor intralayer and interlayer tunneling. It nevertheless suggests the 
likelihood of interesting interaction physics and broken symmetry ground states in 
many neutral or weakly doped multilayer graphene samples. 
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